或许怎么做都对呢？
昨天我们分享了 难者不会，会者不难的心得体会，引起了大家的共鸣，行业内的小伙伴理解了我们的生信服务的正当性，外行人也认可了我们的价值。本来是文末凑数的封面图反而被大家留言纷纷讨论，有小伙伴指出来了这个 Cancer Cell 2022 May文章《Ovarian cancer immunogenicity is governed by a narrow subset of progenitor tissue-resident memory T cells.》里面的差异分析后的，上下调基因，居然是分开多生物学功能数据库注释，然后留言：「最后一个图，图C中，将UP和DOWN 分开分析，是否合理呢？是否合在一起分析更加合理？」
我看了看图例：(C) Gene Ontology (GO) enrichment analysis for differentially expressed genes in TRM-like versus recirculating CD8+ TILs from the same tumors.
比较意外，因为我们一直以来默认授课就是让大家首先差异分析，拿到统计学显著的上下调基因，然后分别注释。我看了看这个数据集：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE194383
是 gene expression profile of CD8+ TILs (TRM vs. re-circulating) from 7 HGSOC ，如下所示样品：
GSM5834878 Darmouth033009_TRM
GSM5834879 Darmouth033009_ReCirc
GSM5834880 Darmouth080609_TRM
GSM5834881 Darmouth080609_ReCirc
GSM5834882 Darmouth092107_TRM
GSM5834883 Darmouth092107_ReCirc
GSM5834884 Darmouth051110_TRM
GSM5834885 Darmouth051110_ReCirc
GSM5834886 042017M_TRM
GSM5834887 042017M_ReCirc
GSM5834888 092817M_TRM
GSM5834889 092817M_ReCirc
GSM5834890 092608M_TRM
GSM5834891 092608M_ReCirc
作者也提供了定量好的表达量矩阵，GSE194383_raw_counts.txt.gz，所以很容易下载并且复现这个过程。但是读者很明显并不是想学数据分析本身，而是讨论分析结果的取舍了。实际上，如果大家使用我的流程，我定义好的 run_kegg 函数 ：
## KEGG pathway analysis
### 做KEGG数据集超几何分布检验分析，重点在结果的可视化及生物学意义的理解。
run_kegg <- function(gene_up,gene_down,pro='test'){
library(ggplot2)
gene_up=unique(gene_up)
gene_down=unique(gene_down)
gene_diff=unique(c(gene_up,gene_down))
### over-representation test
# 下面把3个基因集分开做超几何分布检验
# 首先是上调基因集。
kk.up <- enrichKEGG(gene = gene_up,
organism = 'hsa',
#universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
head(kk.up)[,1:6]
kk=kk.up
dotplot(kk)
kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
write.csv(kk@result,paste0(pro,'_kk.up.csv'))
# 首先是下调基因集。
kk.down <- enrichKEGG(gene = gene_down,
organism = 'hsa',
#universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
head(kk.down)[,1:6]
kk=kk.down
dotplot(kk)
kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
write.csv(kk@result,paste0(pro,'_kk.down.csv'))
# 最后是上下调合并后的基因集。
kk.diff <- enrichKEGG(gene = gene_diff,
organism = 'hsa',
pvalueCutoff = 0.05)
head(kk.diff)[,1:6]
kk=kk.diff
dotplot(kk)
kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
write.csv(kk@result,paste0(pro,'_kk.diff.csv'))
也就是说，我对上下调基因独立做了生物学注释的同时，其实也做了它们合并的，但是我在绘图的时候，仅仅是选取了上下调基因的注释结果，所以会得到文章类似的 图 ：
kegg_plot <- function(up_kegg,down_kegg){
dat=rbind(up_kegg,down_kegg)
colnames(dat)
dat$pvalue = -log10(dat$pvalue)
dat$pvalue=dat$pvalue*dat$group
dat=dat[order(dat$pvalue,decreasing = F),]
g_kegg<- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) +
geom_bar(stat="identity") +
scale_fill_gradient(low="blue",high="red",guide = FALSE) +
scale_x_discrete(name ="Pathway names") +
scale_y_continuous(name ="log10P-value") +
coord_flip() + theme_bw()+theme(plot.title = element_text(hjust = 0.5))+
ggtitle("Pathway Enrichment")
}
kegg_diff_dt <- as.data.frame(kk.diff)
kegg_down_dt <- as.data.frame(kk.down)
kegg_up_dt <- as.data.frame(kk.up)
down_kegg<-kegg_down_dt[kegg_down_dt$pvalue<0.01,];down_kegg$group=-1
up_kegg<-kegg_up_dt[kegg_up_dt$pvalue<0.01,];up_kegg$group=1
#画图设置, 这个图很丑，大家可以自行修改。
g_kegg=kegg_plot(up_kegg,down_kegg)
print(g_kegg)
ggsave(g_kegg,filename = paste0(pro,'_kegg_up_down.png') )
并不意味着 上下调基因合并后的差异基因，它作为一个整体的注释就不重要。我们做数据分析授课的时候，其实更多是方法学的演练，但是结果的解读，更多是生物学背景了。
这样的生物学背景其实绝大部分硕博士以及博士后都不缺的，但是他们如果没有深入了解过前面的上下调基因富集图表的制作流程：表达量矩阵 --> 分组 --> 差异 --> 统计学显著上下调 --> 超几何分布检验 --> 数据库里面统计学显著通路，就没办法把分析结果往生物学意义上面靠。
所以会纠结，到底是将UP和DOWN 分开分析呢，还是合在一起分析更加合理，这样的问题。也会纠结一个通路居然在上下调基因集里面都富集到了，但是，真正理解数据分析过程的小伙伴就会很容易理解图表，并不一定要自己亲自动手分析数据得到图表，熟练一下这个过程（表达量矩阵 --> 分组 --> 差异 --> 统计学显著上下调 --> 超几何分布检验 --> 数据库里面统计学显著通路）就挺好的。
文末友情宣传
强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：
数据挖掘（GEO,TCGA,单细胞）2022年5~6月场，快速了解一些生物信息学应用图表 生信入门课-2022年5~6月场，你的生物信息学第一课